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Compton camera-based prompt gamma (PG) imaging has been proposed for range verification during proton 
therapy. However, a deviation between the PG and dose distributions, as well as the difference between the 
reconstructed PG and exact values, limit the effectiveness of the approach in accurate range monitoring during 
clinical applications. The aim of the study was to realize a PG-based dose reconstruction with a Compton 
camera, thereby further improving the prediction accuracy of in-vivo range verification and providing a novel 
method for beam monitoring during proton therapy. In this paper, we present an approach based on a subset- 
driven origin ensemble with resolution recovery and a double evolutionary algorithm to reconstruct the dose 
depth profile (DDP) from the gamma events obtained by a Cadmium—Zinc—Telluride Compton camera with 
limited position and energy resolution. Simulations of proton pencil beams with clinical particle rate irradiating 
phantoms made of different materials and the CT-based thoracic phantom were used to evaluate the feasibility of 
the proposed method. The results show that for the monoenergetic proton pencil beam irradiating homogeneous- 
material box phantom, the accuracy of the reconstructed DDP was within 0.3 mm for range prediction and within 
5.2% for dose prediction. In particular, for 1.6-Gy irradiation in the therapy simulation of thoracic tumors, the 
range deviation of the reconstructed spread-out Bragg peak was within 0.8 mm, and the relative dose deviation 
in the peak area was less than 7% compared to the exact values. The results demonstrate the potential and 
feasibility of the proposed method in future Compton-based accurate dose reconstruction and range verification 
during proton therapy. 


Keywords: Prompt gamma imaging, dose reconstruction, range verification, origin ensemble, Compton camera, evolutionary 


algorithm. 


I. INTRODUCTION 


Proton therapy has been rapidly developed and widely used 
in clinical cancer treatment over the past decades [1]-[3]. The 
Bragg peak of the proton beam makes it possible to accurately 
deliver enough dose to the target tumors and reduce the radia- 
tion damage to healthy tissues. However, the uncertainties of 
the in-vivo range and dose, which are caused by the treatment 
plan, patient positioning, tumor movement, etc., restrict the 
full clinical potential of proton therapy [4][5]. Range verifica- 
tion is key for further improving the clinical effectiveness of 
proton therapy. The gamma rays derived from proton-induced 
excited nuclei (e.g., !?C* and 160*) are almost prompt emis- 
sion (less than 107+! s); they are called prompt gamma (PG), 
whose distribution is highly correlated with in-vivo dose dis- 
tribution [6]. Using PG for range verification is a feasible 
approach that has been proven in clinical applications [6]. 

A Compton camera (CC) exploits the electronic collima- 
tion principle to realize imaging of radioactive sources. It 
has a higher detection efficiency and abilities of multi-energy 
reconstruction and three-dimensional imaging compared to 
gamma detectors with a mechanical collimator such as single- 
photon emission computed tomography (SPECT). Because of 
its unique advantages, CCs have been proposed for PG imag- 
ing during proton therapy. Several studies on CCs for PG 
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imaging have been conducted, including on Compton sys- 
tem design [7]-[10] and optimization [11]-[13], as well as 
reconstruction algorithm optimization [14]-[18] and novel al- 
gorithms [19]-[22]. Besides, range verification based on CC 
imaging has been experimentally verified [23][24]. However, 
the current range verification based on CCs can only predict 
the peak value and distal falling-off of the PG distribution; 
the prediction accuracy of the actual dose coverage area and 
its relative value is low. If dose reconstruction can be realized 
directly, it will boost the effect of this technology in practi- 
cal clinical applications. However, there are few studies on 
Compton-based dose prediction for proton therapy. The dif- 
ficulty of this problem mainly involves two aspects: 1) the 
accuracy of PG reconstruction is limited by the resolution 
metrics of the detectors and reconstruction algorithm; 2) cal- 
culating the dose distribution from the PG distribution is a 
challenging problem. 


Parodi and Bortfeld proposed a filtering approach to obtain 
the positron emission computed tomography (PET) images 
from dose distribution, assuming that the distribution of the 
positron emitter is the convolution of a pre-defined filter func- 
tion and the dose distribution [25]. The convolution formal- 
ism based on functions was introduced to deduce the dose dis- 
tribution from the secondary radiation distribution. Besides, 
several modified deconvolution methods have been investi- 
gated and evaluated [26]-[28]. Schumann et al proposed a 
novel deconvolution method based on the evolutionary algo- 
rithm (EA) to deduce the dose depth distribution from the PG 
depth profile [29]. Both deconvolution approaches require 
a pre-defined filter function, which depends on the projec- 
tile and target. Another approach is deep learning [30]-[32]. 
Given that deep learning requires a large number of train- 
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ing samples, the results could be unpredictable for unknown 
patients’ bodies and different organizations. By contrast, the 
deconvolution approaches could be more portable. 

The goal of the study was to realize a PG-based dose re- 
construction with a non-ideal CC, thereby further improving 
the prediction accuracy of in-vivo range verification and pro- 
viding a novel method for beam monitoring during proton 
therapy. A modified subset-driven origin ensemble with res- 
olution recovery (SD-OE-RR) is proposed to realize the PG 
reconstruction. We also propose a double evolutionary algo- 
rithm (DEA) to reconstruct the dose depth profile (DDP) from 
the reconstructed PG. We simulated proton pencil beams irra- 
diating into the phantoms made of different materials and the 
CT-based thoracic phantom, respectively. A non-ideal two- 
layer Cadmium—Zinc—Telluride (CZT) CC was used to de- 
tect the proton-induced PG. Finally, the proposed approach 
was evaluated by comparing the reconstructed dose distribu- 
tion with the exact values obtained by the treatment planning 
system (TPS). 


I. METHODS 


A. Monte Carlo Simulation 


Geant4 version 10.03.p01 and GATE version 9.0 with the 
QGSP_BERT_HP_EMY physics list were used for proton- 
induced PG emission and CC detection simulation. Besides, 
the G4EMLivermorePhysics list was used to simulate the 
physics processes in the production and interaction effect of 
PG, including the Doppler broadening effect [33]. The beam 
time structure was referred to the IBA cyclotron C230 used 
in the clinical proton therapy process [34]. The intensity of 
C230 clinical treatment was approximately 2 x 101° s71, in 
which the current was approximately 3.2 nA. The beam pulse 
duration was 3.2 ns with a period of 9.4 ns. In this case, 
the number of protons contained in a single pulse was 217. 
The approximate relationship between the number of protons 
transported N, and the delivered dose is given by Equation 


(1) [23]. 


D 


N, = 6.24 x 10° 
ý S/p 


A, () 


where D is the delivered dose expressed in Gy(J- kg~*) 
and A, is the beam area expressed in cm?. For a 145-MeV 
proton pencil beam with a 2D Gaussian broadening of o=5 
mm, A, was approximately 1.08 cm?. Note also that S/p is 
the mass proton stopping power expressed in MeV: cm?- g~!. 
The average energy of the proton beam at the depth of the 
Bragg peak was approximately 14 MeV, corresponding to a 
residual range of approximately 2 mm. According to the Re- 
port 90 of the International Commission on Radiation Units 
and Measurements (ICRU) [35], S/p was approximately 35.2 
MeV: g71. cm~?. When the delivered dose at the Bragg peak 
was 2 Gy, the number of protons transported calculated by 
Equation (1) was 3.8 x 108. 
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A non-ideal two-stage CZT CC (i.e., presenting limited po- 
sition and energy resolution) with a coincidence time window 
of 1.5 us was used to detect PGs. Each stage of the CZT de- 
tectors comprised 44 x 44 crystals, each with a size of 2 mm 
x 2mm x 15 mm. The spatial resolutions of the detectors 
were both 1 mm in the lateral and depth directions. The en- 
ergy resolution was set as 1.5% at 662 keV. The time resolu- 
tion of the CC was set as 25 ns and the dead time of each event 
was approximately 250 ns. The CC recorded the deposited 
energies and positions of the pixel where the photons inter- 
acted in the two layers and output the time-series projection 
that contained a timestamp and the corresponding interaction 
pixel label and deposited energy. The list-mode data were 
obtained by selecting the time-series events in a coincidence 
time window in which the following requirements were met: 
(i) the total deposited energy was the characteristic energy of 
PGs; (ii) the interactions of the events were once with the 
scatterer first and once with the absorber. The ideal detection 
efficiency of the simulated CC was approximately 1 x 107° 
PG coincidence events for one irradiated proton. However, 
the true coincidence yield, which would determine the de- 
tection efficiency of effective events in practice, was affected 
by the limited coincidence time window of the CC and dose 
rate of the proton beam. The dose rate could be calculated 
by the particle rate of protons using Equation (1). After con- 
sidering the incorrect coincidence caused by multiple protons 
in a bunch and different secondary particles in a coincidence 
time window, the true coincidence yield of the CC for differ- 
ent particle rates of protons was evaluated by simulation; the 
results are provided in Table 1. 


Table 1. Relationship between the particle rate of protons and true 
coincidence yield of the simulated two-layer CZT CC. 


Particle rate of protons particles/bunch true coincidence yield 


2x 10! s~} 217 <0.1% 
2x 10°? s7! 22 28% 
3x107 2 89.5% 
6.7 x 10% s71 1 98% 


Given that the proton beam in a bunch was too large in the 
case of the delivered current in the range of several nA, the 
probability of incorrect coincidence events when using the 
simulated CC detection was greater than 99.9%. Therefore, 
it is almost impossible to obtain correct coincidence events 
that can be used for reconstruction. One feasible method is to 
reduce the current when delivering the beam and prolong the 
irradiation time. When the current intensity was reduced by a 
factor of 10, the number of protons in a single reactor would 
be reduced to approximately 22 by using a delivered current 
of approximately 0.32 nA. For the two-layer CZT CC, the 
true coincidence yield was approximately 28%. After reduc- 
ing the current intensity by two orders of magnitude to 0.03 
nA, the number of protons in a single bunch was 2, and the 
true coincidence yield was approximately 90%. The parti- 
cle rate delivered was 2 x 108 s71, and the time required to 
deliver a 2-Gy dose was 1.9 s. With this particle rate, the rela- 
tionship between the delivered dose, number of protons, and 
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irradiation time used in simulations is listed in Table 2. 


Table 2. Relationship between the total number of delivered protons, 
dose, and irradiation time with the beam time structure used in the 
simulations. 


Number of protons dose irradiation time 
1 x 10° 5.3 cGy 50 ms 
5 x 107 26.5cGy 250 ms 
1x 10° 53 cGy 500 ms 
3 x 108 1.6 Gy 1.55s 
1x 10° 5.3 Gy 5s 


Fig. 1 shows the diagrams resulting from the Monte 
Carlo simulations. To evaluate the performance of the 
SD-OE-RR algorithm for PG distribution reconstructions, 
a 120-MeV proton beam irradiated the water box phantom 
three times independently. Then, a proton pencil beam with 
different energies irradiated the box phantom composed of 
different body-like materials to evaluate the proposed dose 
reconstruction approach. The materials (e.g., muscle, cortical 
bone, soft tissue) used in the simulations were referred to 
the ICRT-46A. To evaluate the performance of the proposed 
method in a situation close to the clinical proton therapy, 
we simulated the proton therapy in the thoracic cavity and 
used the proposed method to reconstruct the dose depth 
profiles with the convolution vectors obtained by the above 
box phantom simulations. The thoracic cavity phantom was 
referred to the Geant4 extended example DICOM [36]. Three 
different proton pencil beams of 2D Gaussian shape in the 
transversal plane irradiated the esophagus or mediastinum in 
different directions to simulate cancer treatment in diverse 
situations. For the thoracic tumor proton therapy simulation, 
two treatment modalities were considered. One was a single 
spot scan of the pencil beam with approximately 0.8-MeV 
energy broadening and spatial broadening parameters o 
equal to 2 mm; the coverage depth of the Bragg peak maxi- 
mum irradiation dose was approximately 1 mm. The other 
modality was the spread-out Bragg peak (SOBP) irradiation 
with 5.8-MeV energy broadening and spatial broadening pa- 
rameter o=5 mm for full coverage of the target area. A total 
of 1 x 107 protons were used in the simulations to evaluate 
the proposed method. Finally, to investigate its performance 
with a different dose, various numbers of protons, ranging 
from 107 to 10°, were implemented for the thoracic cavity 
proton therapy. The corresponding doses of the different 
numbers of protons delivered to the phantoms varied from 


oij = In[{1 + d6(max{| cos 64; — cos 0;| — A(cos 44;)|, 0})]/In2 


f (P1041) = max{sen(2(F(61,441) + 1— F(G1,x))), sen(2(F' (61,4) — F(81,%+1))) > sen( 


- (F(61,441) + 1], sgn(2(F (81,4 — F(G1,441)))) - sgn( 2,441 — 
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5.3 cGy to 5.3 Gy, covering the region of generally delivered 
dose in clinical treatments [37]. 


B. Subset-driven origin ensemble with resolution recovery 


The subset-driven origin ensemble with resolution recov- 
ery (SD-OE-RR) algorithm was used for PG reconstruction 
with list-mode projection data. Different from the SD-OE- 
RR proposed in a previous study [38], the calculations of the 
resolution correction factor and initial guess of the source dis- 
tribution fo given by Equations (2) and (3) were modified for 
CZT CC-based PG reconstruction. 


2 
Moc 
A(cos bgi) S in ‘(nei + NDi) + NPI 


Ta (2) 


fo = 45, 4(v;, ei) | cos Ori — cos 0;| < A(cos 0gi)} (3) 


where cos @; denotes the cosine of the Compton scatter- 
ing angle calculated by the deposited energies for event 2; 
cos fs is the theoretical cosine of scattering angle determined 
by the scattering position, absorbed position, and voxel v; in 
the field of view (FOV); E;2 is the deposited energy in the 
absorber; np; is the deviation caused by the spatial resolution 
of the detectors; and 7”; and 7p; are the energy deviation 
percentages caused by the statistical fluctuation and Doppler 
broadening effect, respectively. For the CZT CC used in the 
study, the total cosine A(cos 4,;) was approximately equal to 
0.05. 

The voxels probably containing the sources of event 7 were 
stored as a one-dimensional vector O; = %,j0;;, where the 
subset of the origin ensemble };0;; is given by Equation (4). 
Then, the iteration of the SD-OE-RR was based on the Monte 
Carlo-Markov chain, updating the probability density func- 
tion by Equation (5); this is similar to the original origin en- 
semble algorithm [19]. The defined (x) function equals 0 
when x Æ 0 and equals 1 when x = 0. Moreover, k+1 and k 
represent iteration times, 31,41 represents the randomly se- 
lected subset of the origin ensemble in the (k + 1)th iteration, 
and (2,441 represents another independent random number 
used in the (k + 1)th iteration to determine whether to move 
the location of the origin. Besides, F (61, %+1) is the sum of 
f(G1,n41) and f of the 81 ,41-centered four adjacent pixels. 
Finally, sgn(x) denotes the sign function. 


(4) 


F(Bi k41) +1 
F(61,k) 


F(61,441) +1 
Fin? [F(81,4%) + 1]} 


— b2,k+1) 
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Fig. 1. Diagrams resulting from the Monte Carlo simulations; (a) CC-based beam range and dose monitoring in the multi-layer box phantom; 
(b) non-ideal CC with limited position and energy resolution; (c) proton therapy simulation of the thoracic cavity based on the CT slice. 


SD-OE-RR was programmed with CUDA C++ for paral- 
lel acceleration. The FOV size was set as 200 mm and 107 
iterations were implemented for each reconstruction. 


C. Double evolutionary algorithm 


In a previous study, the PG depth profiles (GDPs) and dose 
depth profiles (DDPs) could be fitted by an analytical approx- 


imation of the Bragg curve called Q function. The convolu- 


tion relation between the corresponding fitting curve GDP 


and DDP was given by Equation (6), which is an analytical 
approximation of the Bragg curve. The Q function is the con- 
volution of a Gaussian function G(x) and a power-law func- 
tion P, (x), which was introduced by Parodi and Bortfeld [20] 
for PET and adapted to PG by Schumann et al [29]. In this 
study, the evolutionary algorithm was used to obtain the ker- 
nel function k of various known simple phantoms and deduce 


the unknown DDP of complex phantoms. This approach was 
abbreviated as double evolutionary algorithm (DEA). 


GDP = DDP -k (6) 
The DEA iteration process used in this study is described 

next: ys 

0. Initialization. A specific array of a Q distribution is used 

to create a set of Npop individuals. Each individual represents 


a DDP ork array for two different evolutionary algorithm 
applications. 

1. Parent selection. Two parental arrays are obtained with a 
fitness proportional selection. The arrays with higher fitness 
give rise to new offspring. 

2. Crossover. With random points cutting in two arrays, the 
parental arrays implement the single-point crossover to form 
two-child arrays. 

3. Mutation. The new offspring is mutated to derive a 
potential improvement in fitness. Each modified array could 
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be shifted in depth by a random integer value in the interval 
[—2,2] with probability p = 0.3 and multiplied with a 
uniformly distributed random factor between [l-r,1+r] 
(r=0.2) with probability pọ = 0.3. It could be also varied 
around a randomly chosen point with probability p = 0.4 
and a Gaussian curve of random height (between +10%) and 
width of 10 mm (20). 

4. Replacement. The next generation is created by choosing 
50% of individuals with the best fitness value from the 
parental and newly generated population. 

5. Iteration. Repeat steps 1 to 4 until reaching the pre- 
determined iterations, 1.e., Niter. 


The specific array mentioned above, which replaces their 
corresponding continuous function, represents a vector with 
400 elements in the interval [-200 mm, 200 mn] for k and a 
vector with 200 elements in the interval [0 mm, 200 mm] for 


DDP, respectively. The DEA iterations were implemented 


3000 times for the kernel function k of various known sim- 
ple phantoms and 1000 times for unknown DDP of complex 
phantoms. 


D. Dose depth profile reconstruction framework 


As shown in Fig. 2, the CC-based dose depth reconstruc- 


tion framework mainly comprises two parts: the GDP ob- 
tained by the CC and the convolution kernel k obtained by 
prior known phantoms and final dose estimation via the DEA. 
Then, the estimated convolution kernel k of multiple materi- 
als or multi-energy was obtained by algebraic averaging (1.e., 
interpolation method) or weighted average of the known ker- 
nels. The previously known kernels with proton beams of 
different energies irradiating different materials were calcu- 
lated by the DEA. The original GDP was obtained by SD- 
OE-RR reconstruction with the projection data of a non-ideal 


CC. Moreover, the estimated GDP was given by local peak 
fitting based on the original GDP. The local fitting interval of 
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Fig. 2. Framework of dose depth profile reconstruction from CC detection-based PG by the proposed method. 


curves used in the study was the interval above 80% of the 
peak value. 


Il. RESULTS 


A. PG reconstruction 


Table 3 shows the PG reconstruction results for three origin 
ensemble (OE) algorithms (i.e., ordered OE-RR [15], ROE- 
RR [18], and SD-OE-RR algorithms). The PGs with four 
characteristic photons (i.e., from 12C, 14N, 150, and 160 
de-excitations) were chosen to evaluate their performance. 
To alleviate the effect due to the incomplete absorption and 
background radiation, the effective events for reconstructions 
were selected by using the total energy windows of coinci- 
dence events within +0.2 MeV of the four known PG energy 
spectral peaks (i.e., 4.44 MeV, 2.31 MeV, 5.25 MeV, and 6.13 
MeV) [15]. As shown in Table 3, compared with previous OE 
algorithms, SD-OE-RR presents the same or slightly higher 
reconstruction accuracy. Besides, in a 64-bit Linux computer 
with a 2.50 GHz Intel i5-7200U CPU and a GTX 1650 Ti 
Nvidia GPU, for approximately 200000 events, the recon- 
struction times were 26, 21, and 3 seconds for the ordered 
OE-RR, ROE-RR, and SD-OE-RR algorithms, respectively. 
Therefore, the SD-OE-RR algorithm proposed in this study 
consumes less time while providing reconstruction whose 
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distal-fall position deviation was less than 2 mm. The results 
also demonstrate that the proposed algorithm can make use 
of the PG distribution reconstructed by the SD-OE-RR algo- 


rithm to obtain the estimated GDP with good agreement in 
the positions of peak and distal fall-off. 


B. Reconstruction of dose depth profiles 


As shown in Figs. 3(a) and (b), compared to the exact DDP 
in the two-layer phantom made of cortical bone and muscle, 
the reconstructed DDP with CC by the proposed method had 
less than 0.3 mm deviation at the 50% distal fall-off position. 
Moreover, the relative deviation was within 2.5% in terms of 
absolute dose in the region above 80% maximum intensity. 
Besides, the reconstructed DDP had an average relative devi- 
ation of 7.7% in the region of the level area behind the dose 
mutation interface. As shown in Figs. 3(c) and (d), the re- 
constructed DDP of the three-layer phantom made of cortical 
bone, muscle, and soft tissue had a deviation below 0.26 mm 
at 50% distal fall-off position, and within 5.2% in terms of 
absolute dose around the Bragg peak. As shown in Figs. 3(e) 
and (f), for the SOBP induced by a 145~160 MeV proton 
beam irradiating a phantom made of cortical bone, muscle, 
and soft tissue, the reconstructed DDP could reproduce the 
peak region with an accuracy within 2 mm and predict the 
relative dose in the peak area with a deviation less than 4%. 


Table 3. Simulated proton-induced PG reconstruction results of the three OE-RR algorithms for 10’ incident protons. 
Position Method  “C “-N OO. ë "o+ o+ Oo 


100.72 + 0.03 


100.86 + 0.03 


Position Method RO 14N 
MC 98.18 + 0.37 85.16 + 0.64 
Peak OE-RR 95.98 + 0.63 81.29 + 0.83 
ROE-RR 97.52 + 0.45 82.53 + 0.87 
SD-OE-RR 96.87 + 0.44 83.52 + 0.93 
MC 90.87 + 0.09 
80% OE-RR 98.12 + 0.68 88.14 + 0.34 
fall-off ROE-RR 99.03 + 0.81 88.34 + 0.53 
SD-OE-RR 99.32 + 0.43 88.96 + 0.88 
MC 94.22 + 0.22 
50% OE-RR 101.97 + 0.63 92.33 + 0.26 
fall-off ROE-RR 101.06 + 0.35 91.95 + 0.92 
SD-OE-RR 101.67 + 0.78 91.85 + 0.52 


BO TO 
95.31 0.39 101.30 + 0.37 100.00 + 0.39 
89.64 + 0.53 100.20 + 0.42 97.95 + 1.07 
93.32 + 0.61 99.96 + 0.37 99.14 + 0.45 
93.67 + 0.56 100.72 + 0.35 99.63 + 0.39 
97.68 + 0.10 103.76 + 0.06 102.75 + 0.09 
93.47 + 0.51 104.12 + 0.15 102.57 + 0.60 
96.31 + 0.46 103.10 + 0.55 102.21 + 0.27 
96.12 + 0.52 103.19 + 0.06 102.26 + 0.36 
98.95 + 0.51 104.06 + 0.06 104.49 + 0.22 
101.39 + 0.34 105.75 + 0.32 104.53 + 0.43 
100.38 + 0.73 105.87 + 0.31 104.85 + 0.73 
99.73 + 0.49 105.57 + 0.25 104.57 + 0.30 


345 Furthermore, the relative deviation of the global curves (until 350 
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the distal end fall-off) was approximately 4.5%. 
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Fig. 3. Simulation results of proton beam irradiation for multi-layer 
box phantoms: (a)(c)(e) show the comparison of the DDPs between 
the exact values obtained by MC simulations and reconstructed val- 
ues via the proposed method from CC for three different multi-layer 
phantoms consisting of cortical bone, muscle, and soft tissue, re- 
spectively; (b)(d)(f) show their corresponding relative dose devia- 
tion, respectively. 


Fig. 4. shows that, with respect to the exact 2D dose dis- 
tributions, the Compton-based reconstructed PG distributions 
obtained by the SD-OE-RR algorithm were in good agree- 
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ment with the positions of the distal fall-off, but the peak 
broadening was skewed, similar to the distributions of the ini- 
tial PGs. Note also in Figs. 4(d) and (f) that the reconstructed 
DDPs for the in-vivo proton beam had a deviation of less than 
0.6 mm for the distal fall-off position. Moreover, the recon- 
structed DDPs were in good agreement with the exact values 
at the Bragg peak in the region above 80% maximum inten- 
sity, where the relative deviation was within 5.3%. However, 
as shown in Figs. 4(d) and (e), the reconstruction in the region 
before the Bragg peak had a large deviation due to the pro- 
ton beam passing through a more complex structure such as 
a lung. In contrast, Fig. 4s(i) and (j) show that when the pro- 
ton beam passes through a more homogeneous structure such 
as mediastinum, most of the relative deviations were within 
10% in the region of the level area behind the dose mutation 
interface. 


Compared to the exact initial PG distribution shown in Fig. 
5(b), the Compton-based reconstructed PG distributions ob- 
tained by the SD-OE-RR algorithm shown in Fig. 5(c) were 
in good agreement with the distal fall-off position and spa- 
tial distribution in the area around the peak. As shown in 
Fig. 5(d), the reconstructed DDP for the in-vivo proton beam 
could predict the position of the distal fall-off with an accu- 
racy within 0.8 mm. Moreover, as shown in Fig. 5(e), the re- 
constructed DDP was in good agreement with the exact value 
at the Bragg peak in the region above 80% maximum inten- 
sity, with a relative deviation of less than 4.8% in terms of 
absolute dose. However, for the region of the level area be- 
fore the Bragg peak, the reconstructed DDP had larger values 
than the exact values with a mean relative deviation of ap- 
proximately 9.2%. 


Fig. 6 shows the simulation results of proton therapy for 
thoracic mediastinal tumors with different numbers of inci- 
dent protons. Given that the reconstruction accuracy of the 
PGs with the CC was influenced by the number of detected 
events, the DDP reconstructed by the proposed method was 
correlated with the number of incident protons. When the 
number of incident protons increased, the number of effective 
events available for the reconstruction of PGs also increased 
by almost the same proportion, under the premise that the par- 
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Fig. 4. Simulation results of proton therapy for thoracic mediasti- 
nal tumors with 130.2~131 MeV proton pencil beams ((a)~(e)) and 
for thoracic esophageal tumors with 144.2~145 MeV proton pencil 
beams ((f)~G)), respectively. (107 protons). (a)(f) show the exact 
distributions of dose at the CT slice obtained from MC simulations; 
(b)(g) show the distributions of the initial PGs induced by the pro- 
ton beams; (c)(h) show the reconstructed PG distributions obtained 
by SD-OE-RR with CC data; (d)(i) show the comparison between 
the exact DDPs in (a)(f) and reconstructed DDPs from (c)(h); (e)(j) 
show their corresponding relative deviations. 


ticle rate of incident protons (which was proportional to the 
dose rate) remained unchanged. When the number of inci- 
dent protons varied from 10” to 10°, corresponding to a deliv- 
ered dose ranging from 5.3 cGy to 5.3 Gy, the number of PG 
events by the simulated CZT CC ranged from approximately 
2000 to approximately 200000. As shown in Fig. 6(a), the 
accuracy of the reconstructed GDP deteriorated as the parti- 
cle number decreased, especially in the area before the Bragg 
peak. This is consistent with the results of previous studies 
on OE-RR-based PG reconstruction [15][18]. However, the 
reconstructed positions of the distal fall-off were almost the 
same, which means that the distal fall-off was reproduced ac- 
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Fig. 5. Simulation results of proton therapy for thoracic mediastinal 
tumors with 131.2~137 MeV proton pencil beams (107 protons) of 
2D Gaussian shape in the transversal plane (o = 5mm): (a) exact 
distribution of dose obtained by MC simulation; (b) corresponding 
initial PG distribution; (c) reconstructed PG distribution obtained by 
SD-OE-RR with CC data; (d) comparison between the exact DDP in 
(a) and the reconstructed DDP from (c); (e) corresponding relative 
deviations of dose. 
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Fig. 6. Simulation results of proton therapy for thoracic mediastinal 
tumors with different numbers of incident protons: (a) reconstructed 
GDPs obtained by SD-OE-RR for different particle rates; (b) relative 
deviations between the exact and reconstructed DDPs for different 
particle rates; (c) corresponding reconstructed DDPs. 


4o2 Curately. Therefore, the reconstructed DDP could provide an 
s3 accurate range prediction for 10” incident protons. Moreover, 
44 the accuracy of the dose prediction will be better if the num- 
45 ber of protons increases. As shown in Figs. 6(b) and (c), 
as for a total 1.6-Gy dose, which corresponds to 3 x 108 in- 
47 cident protons, the reconstructed DDPs were in good agree- 
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ment with the distal falling edge of the Bragg peaks with an 
accuracy within 0.6 mm, and the corresponding relative de- 
viations in the region above 80% maximum intensity were 
less than 5.5%. When the total delivered protons reached 10° 
(i.e., approximately 5.3 Gy), the relative deviations in the re- 
gion above 80% maximum intensity were less than 4% and 
the average value of the global deviation was approximately 
5%. 


IV. DISCUSSION 


The difficulty of dose reconstruction based on CC comes 
from two aspects: the rapid and accurate PG reconstruction 
with non-ideal CC measured data and the calculation of dose 
from the reconstructed PG distribution. In this paper, we first 
propose a modified SD-OE-RR algorithm to realize a faster 
and more accurate PG reconstruction. As shown in Table 3, 
the SD-OE-RR algorithms provided peak estimations with an 
accuracy of approximately 1.0 mm for all the PGs. The accu- 
racy of the distal fall-off positions was within 2 mm except for 
the 50% fall-off position of £N. Compared with the previous 
OE-RR algorithm, the proposed SD-OE-RR algorithm addi- 
tionally considered the correction of the Doppler broadening 
effect calculated by the CZT extranuclear electron momen- 
tum, which further improves the speed of the convergence and 
accuracy of the relative intensity peak of the reconstructed PG 
[39]. Besides, the parallel SD-OE-RR algorithm exploits the 
advantages of GPU multi-thread simultaneous computing and 
greatly reduces the reconstruction time, meeting the require- 
ments of rapid reconstruction with a large number of events. 
For the proposed SD-OE-RR algorithm, in order to reduce 
the image blur caused by the large variability between the av- 
erage states for computing the reconstructed images during 
OE iterations [17], we selected the imaging field of view with 
a number of voxels of 256 at most in each dimension while 
ensuring that it covered at least 20 cm in the depth direction 
to achieve beam monitoring. Besides, we only used the total 
energy window to select the effective events of PGs. How- 
ever, we identified sources of background noise, such as sec- 
ondary protons or neutrons generated during the irradiation, 
that could lead to an accidental coincidence event of wrong 
origin in the actual irradiation; this could hardly be culled 
by the energy windows [40]. Better event selection methods 
such as those based on neural networks [41] or the distance- 
of-closest approach [42] can be used for PG event selection 
before reconstruction by the proposed method in future clin- 
ical applications. After PG reconstruction, the corresponding 
profile (i.e., GDP) with SD-OE-RR for 160 was used to ob- 


tain the GDP by curve local fitting with Q. Selecting the 
fitting interval above 80% of the peak value was due to the 
good agreement between the reconstructed and exact values, 
as shown in a previous study for OE-RR-based PG recon- 
struction. a 
Based on the original evolutionary algorithm and Q func- 
tion fitting proposed in Parodi and Bortfled [25], we pro- 
pose a DEA algorithm to realize the dose calculation from 
the reconstructed PG distribution. The main modifications 
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were the calculation approach of the convolution kernel and 
a novel convolution kernel calculation approach for complex 
phantoms. Instead of using the analytical method, we used 


a similar EA to obtain the vectors k with the same length 
of elements to replace the continuous kernel function. We 
implemented the protons with an energy range from 131 to 
159 MeV while irradiating different phantoms. In this case, 
the main peak positions of these vectors remained almost un- 
changed when the energy of the proton beams or the density 
of targets did not vary significantly. Given that the SOBP gen- 
erally uses an energy range within 30 MeV, the filter kernel of 
the corresponding proton energy in the distal falling edge or 
median depth can be used for dose prediction, as verified by 
Schumann et al. [29]. The density of most human tissues is 
between 1.0~1.1 g/cm, except for the cortical bone, which 
is approximately 1.82 g/cm. Thus, the average convolution 
kernel could be calculated from the convolution vectors ob- 
tained by the monoenergetic proton beams irradiating the box 
phantoms composed of a single material. Moreover, the cor- 


responding estimated convolution vector k were used as the 
pre-defined filter kernel in multi-layer materials and polyen- 
ergetic protons, as well as the complex thoracic phantom. In 
particular, because the specific material proportion of the tho- 
racic phantom and their convolution kernel vectors were often 


not known in advance, its average convolution vector k was 
calculated using the convolution kernel vectors of materials 
with similar densities (i.e., water and cortical bone), to simu- 
late the practical application of the proposed method. 

The results shown in Fig. 3~Fig. 6 verified the feasibility 
of the proposed method in the range and dose prediction dur- 
ing proton therapy. The DDPs reconstructed by PGs with a 
non-ideal CC were in good agreement with the exact values of 
dose distribution calculated by TPS. In the simulation of the 
box phantom, the reconstructed DDP had an accuracy of less 
than 0.3 mm for range prediction and within 5.2% for dose 
prediction of monoenergetic protons irradiating. The predic- 
tion deviation of the range was less than 2 mm for a polyener- 
getic proton beam irradiating the multi-layer phantom (even 
for a large spread-out Bragg peak of approximately 2 cm). 
Besides, for the proton pencil beam spot scanning simulations 
with SOBP of approximately 3 mm with a total dose (approx- 
imately 2 Gy ~5 Gy) commonly used in clinical applications, 
the reconstructed DDPs for the in-vivo proton beam had less 
than 0.8 mm deviation for the distal fall-off position and were 
in good agreement with the exact values in the region above 
80% maximum intensity. The accuracy of the reconstruction 
depends on the number of incident protons, which determined 
the statistical properties of the measured projection data. The 
deviations of the range prediction were less than 1 mm for a 
number of incident protons of at least 10’, which corresponds 
to a delivered dose of 5.3 cGy. Moreover, the accuracy of the 
dose prediction is better as the number of incident protons 
increases. 

Compared with the previous range prediction methods 
based on CC-based PG imaging, the proposed method 
achieved higher accuracy in terms of range verification by 
directly reconstructing the dose depth distribution. This is 
because the distal falling edge of the PGs and that of the dose 
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distribution usually had a deviation of approximately 2 mm 
or more [43]. Therefore, the indirect range prediction based 
on the PGs had systematic errors. The proposed method si- 
multaneously realized the dose reconstruction of the Bragg 
peak region of the proton beam in different materials, thereby 
providing a means for online monitoring and control of the 
beam hot spot. Compared with other offline dose verifica- 
tion methods, the proposed dose prediction method based on 
CC is promising in terms of optimization of the parameters of 
the proton beam and therapy plan according to the patient’s 
condition after a short-time feedback within several seconds, 
and to further improve the curative effectiveness of proton 
therapy. In future work, after completing the construction of 
the CC prototype, we will conduct range and dose prediction 
experiments based on CC applied on the proton beam under 
clinical conditions to further optimize the proposed method 
and apply it in future rapid beam monitoring. 


V. CONCLUSIONS 


In this paper, we propose an approach to realize the dose 
depth profile reconstruction with a CC for proton therapy. A 
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modified SD-OE-RR algorithm and double evolutionary al- 
gorithm are presented to address the two main challenges in 
CC-based dose reconstruction. The simulation results of the 
box phantom irradiation and thoracic tumor treatment demon- 
strate the feasibility of the proposed method in accurate dose 
depth distribution reconstruction. Moreover, the simulation 
results also demonstrated that the algorithms proposed in this 
paper are feasible for more accurate prediction of beam range 
based on CC during proton therapy. The proposed method 
may be used in future rapid dose monitoring and more accu- 
rate range verification during proton therapy. 
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